Shock waves in ultracold Fermi (Tonks) gases 
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It is shown that a broad density perturbation in a Fermi (Tonks) cloud takes a shock wave form in 
the course of time evolution. A very accurate analytical description of shock formation is provided. 
A simple experimental setup for the observation of shocks is discussed. 



Shock waves show up in various physical systems. For 
example, we find them in a gas bubble driven acoustically, 
in a photonic crystal, and even in a Bose-Einstein con- 
densate Q. In this Letter we study properties of shock 
waves in ultracold Fermi (Tonks) gases: the simplest 
many-body quantum system. To this aim, we consider 
dynamics of a Gaussian-like density perturbation being 
initially at rest on the fermionic cloud. We show that 
such a perturbation divides into two pieces that travel in 
opposite directions in the course of free time evolution. 
The two perturbations change their shape during propa- 
ation, and eventually two shock wave fronts are formed 
We provide a detailed theoretical description of their 
formation and further propagation. 

Experimental studies of such phenomena are possi- 
ble due to recent progress in cooling and trapping of 
fermionic gases 0|. What's more, it turns out that inter- 
particle interactions between cold trapped fermions are 
neglig ible 0, making the theoretical description exactly 
tractable. 

Another exactly solvable many-body system is related 
to one-dimensional (ID) hard-core bosonic gas, the so 
called Tonks gas 0, la, @ • As suggested by Olshanii et 
al. || tightly confined alkali atoms are very promising 
candidates for experimental realization of the Tonks gas. 

There is a close connection between Tonks and ID 
Fermi gases due to the Fermi-Bose mapping theorem 
(FBMT) 0,0. This theorem says that if #(xi, • • • , x N ) 
is the ground state of the A^-fermion system described by 
(dimensionless) Hamiltonian 
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then the exact many-body ground state Tonks wave func- 
tion, in an external potential the same as in is 
ni<j<fc<iv sign(ar A - Xj)*(xi, • • • ,x N ) 0,0. As a re- 
sult, the single-particle density 



Qn{x) = / dx 2 ■ ■ ■dxjv|*(x,x 2 , • • • ,Xjv)| 2 



(2) 



is identical for Fermi and Tonks clouds. Since FBMT 
can be extended to time dependent problems [|| , density 
perturbations on Fermi and Tonks clouds propagate in 
exactly the same way. 

It was proposed by Kolomcisky and Straley that quan- 
titative features of Tonks (Fermi) density profile can be 



correctly addressed by the mean- field approach Q , where 
density qn{x) ~ |<f>(x)| 2 with $(x) being the ground 
state solution of 

" \^ {x) +V ^ X ) + § l*^)! 4 *^) = ( 3 ) 

with G = tt 2 N 2 , N 3> 1, and ^(x)! 2 normalized to unity. 
On the basis of © we find a distance over which the 
density profile tends to its bulk value p when subjected to 
a localized perturbation. This characteristic length scale 
is called the healing length £(p), and equals \j\fGp. A 
natural extension of © allows for time dependent studies 
according to 
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Ydx 2 



V+^)$(x,t). (4) 



Therefore, a many-body character of the problem has 
been overcome by introducing nonlinearity into the the- 
oretical description. 

In this Letter, we investigate formation of shocks and 
their dynamics using both the mean-field approach Q 
and the exact many-body treatment 0. We show that 
the mean-field formalism allows for accurate analytical 
determination of all the quantities concerning shock for- 
mation. Explicit comparison of the mean-field results to 
the exact calculations allows us to specify the range of 
applicability of Eq. <0J. After shock creation the mean- 
field approach becomes inaccurate, and our discussion is 
based solely on exact findings. 

To begin, the dynamics of an initial density perturba- 
tion in a ID system will be described. This perturbation 
can be produced by an adiabatic focusing of an appropri- 
ately detuned laser beam on the particles 0- Alterna- 
tively, one may cool the cloud in the presence of a laser 
beam . We assume that our system consists of a zero 
temperature Fermi (Tonks) gas. Atoms are placed in a 
ID periodic box having boundaries at x = ±Z. 

In the Thomas-Fermi regime, when the kinetic energy 
term in is negligible, density profile takes form 



|$(x,0)| 2 = poV / l + 2ae-* 2 /2<T 2 j 



(5) 



where we assume that the laser produces potential 
Mo exp(— x 2 /2cr 2 ), with £(po) < < i, and a cx «o- 
We determine po from the condition J dx|3>(x, 0)| 2 = 1. 
For a < 0.5 density profile can be rewritten as follows 



|$(x, 0)| 2 w po + po(VT+2c^ - l)e- x2 t 2K * 



(6) 
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where k = a J dx(Vl + 2ae- x2 - l)/(y/l + 2a - 1)Vtt 
provides a correct normalization of (jBJ- 

Suppose that the laser beam has been abruptly turned 
off and the perturbation © starts to evolve. If the sys- 
tem was described by a single-particle Schrodinger equa- 
tion, the perturbation would spread out. To find out 
what happens in a many-body problem first we consider 
the low amplitude limit. The "wave function" is written 
as ®(x,t) = e-^*[$ + 8$(x,t)], where <J> = 1/V21 and 
/i = N 2 ir 2 /8l 2 . The density profile reads 



Mx,t)\ 2 = <f>2 + $ [«5<i>(x,i) + ^(x,t)] + O(|^| 2 ). (7) 
For a-»0we find, with the help of JSJ and Q, that 

6$(x, 0) = a<S>o(e- x2/2a2 - aV2^/2l)/2 + 0{a 2 ). (8) 
The linearization of Q leads to the equation 
1 d 2 



2dx 2 



5$ + G$o(^ + 5®*), 



which can be solved by means of Bogoliubov substitution: 
6$(x,t) =Y, n a n K(x)e- iUnt +Y* n (x)e iu ^}. Straightfor- 
ward calculation determines modes u„, v„ and frequen- 
cies u n . Coefficients {a n } are found from the initial con- 
dition (JSJ. Substitution of 5$>(x, t) into (7J gives 



\$(x,t)\ 2 = ^ 2 Q + aa^ 2 V2^ 



^0 

E 

71=1 



(9) 



e -fc> 2 /2 cos (fc n2 ;) cos(uj n t)/l + 0{a 2 ), 



where ui n = k n y/c^ + k%/2, c = VG® 2 ,, and k n = nn/l. 
Approximating oj n w fc„co (cr ^> 1/co), and changing 
summation into integration we arrive at 



Mx,t)\ 2 



% - a%V2n<j/2l 



a$ 2 [e 



+ e 



-(a:+cot) 2 /2CT : 



]/2 + 0(a 2 ). 



-(x-c t) 2 /2a 2 

(10) 



Thus, we see that the initial density profile breaks into 
two separate pieces moving with the sound velocity Cq 
in the opposite directions. The dynamics of these per- 
turbations will the subject of subsequent considerations. 
Finally, notice that the whole discussion applies to sys- 
tems with N > 1 (see ©). 

Having in mind the low amplitude result, a hydrody- 
namical approach, that fully accounts for nonlinear char- 
acter of the problem, will be developed. We write "wave 
function" $ as y / pexp(ix) and define the velocity field: 
v = d x x- Eq. j3J in (/?, v) variables for a homogeneous 
system takes form 
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FIG. 1: Density of atoms as a function of time. Dotted line 
corresponds to Pn(%, 0). Solid line is the exact JV-fermion 
solution [3 , and a dashed line is the analytical approximation 
itn . Plot (a) [(b)] is for t = 22.2 [44.4]. Parameters: a = 
0.12, Wo ~ -1-16, o = 12.5, k » 12.7, N = 399, I = 200, see 
fll| for units. 



where (|ll|l is a continuity equation, while (|12|) is similar 
to the Euler equation from classical hydrodynamics. To 
simplify the problem we will neglect in l|12fl the so called 
quantum pressure (QP) term ~ d 2 \fp/\fp~ |9j- Simple 
estimation shows that it is of the order of 1/s 2 with s 
being the typical spatial scale of variations of the density 
profile. Therefore, it is important when p 2 G ~ 1/s 2 , 
which means that s is of the order of the healing length £. 
We assume that initially a perturbation is broad: s £. 

Solution of a similar problem in classical hydrodynam- 
ics is known for years Q. Following [3] we find exact 
traveling wave solutions of 1|11|) and l|12l) (without the 
QP term) 



p = f{x - (a ± 2VGp) t) 



no 



where / (an arbitrary function) and ao (a constant) can 
be found from initial conditions. The sign ± determines 
a direction of propagation. This solution is qualitatively 
the same as the one obtained for the ideal ID gas @], 
where the existence of shocks is well known. 

We aim at the determination of dynamics of a single 
Gaussian-like profile moving to the right. Analytical so- 
lution of a full problem seems to be very difficult due to 
the lack of superposition principle. On the basis of the 
low amplitude result H10|) . we assume that initially the 
right moving part is a Gaussian of half-amplitude of the 
initial perturbation ©, and impose a natural condition 
that v(x -> ±1, 0) -> 0. It yields to 



p(x,t) = po+p J]eM-[x-VG(~p +2p)tY /2k 2 ), (13) 



with r] = (VI + 2a - l)/2 ©. We would like to stress 
that l|13|) holds for an arbitrary large r\ since we fully 
account for the nonlinear term. Many interesting results 
can be obtained from i|13|) despite its implicit form. 
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First, the maximum of the impulse has a constant am- 
plitude equal to po(l + rj) and moves with a constant 
velocity V(rj) = vGpo(l + ^v)- F° r 1 ^ we get 
Vo ~ VGpo = Co, which is the speed of sound propa- 
gating at the background density po- I n this limit we 
recover (fTOfl . 

Second, bright perturbations {rj > 0) move with higher 
velocity, V(rj), than dark structures (rj < 0). This sim- 
ple prediction could be directly verified experimentally. 
Moreover, it is in a qualitative agreement with numerical 
simulations of Karpiuk et al. \l2l , even though they 
generate excitations in a quite different way. 

Third, the shape of the perturbation changes in the 
course of time evolution - Fig. ^ Since the speed of 
propagation increases with density, the upper part of 
bright impulse moves faster than its tail leading to self- 
steepening in the direction of propagation. Eventually, 
shock will be formed, i.e. dp{x)/dx = — oo in at least 
one point of the profile. 

To determine time and position of shock formation we 
solve the set of equations [2| 



dp 







d 2 x(p) 
dp 2 



0. 



(14) 



The first one is equivalent to \dp(x)/dx\ — oo, and the 
second one is necessary to assure uniqueness of p(x). For 
the density profile (|T3Tl one easily finds x(p) and solves 
()14|l . Without the loss of generality we assume from now 
on that rj > 0, and find time t s of shock formation, and 
density p s at which shock appears 



ts 



2^/Gp^e- 1 / 2 



Ps = Po + PoV e 1/2 • (15) 



As a result, an initially broad Gaussian- like density per- 
turbation moving in the Fermi (Tonks) cloud (|13(l devel- 
ops a shock wave front during time evolution. This is 
the central result of our paper. The solution (|10[1 fails 
to predict such behavior, since it is derived in the zero 
amplitude limit {rj — > 0) where t s — ► oo. 

Interestingly, the time of shock formation can be cor- 
rectly estimated without solving 1|14(1. Indeed, the half- 
width of the impulse (w 2k) is a difference in a dis- 
tance traveled by the upper and lower parts of the im- 
pulse until shock appears: 2k w (V(ri) — Vo)t s . It gives 
t s k/ \/Gr]po- Estimating the position of a shock wave 
front by position of impulse's maximum we get a very 
simple expression: x s — Ke}l 2 (\ + \j2r\). 

Comparison of the analytical solution l|13Jl to the TV- 
body one || for a moderately large initial perturbation 
(x s > I) shows excellent agreement-see Fig. Q] As x s be- 
comes smaller than half-box-size (= I) shock waves show 
up in our simulations. We have increased the amplitude 
of the initial perturbation and compared a full mean- field 
calculation with the exact fermionic one ||. As shown 
in Fig. |3 until the expected moment of shock wave cre- 
ation agreement is very good. The analytical solution 




FIG. 2: Density of atoms as a function of time. Subsequent 
profiles correspond to t = 41, 82, 131.2. Dashed line is an 
exact TV-body calculation § , while solid line presents numer- 
ical solution of mean- field equation Q ■ Inset shows details of 
a density profile at t — 82. Shock's characteristics: t s = 81.6 
and x s = 118.4. Parameters: a = 0.5, uq ~ —0.53, a — 20, 
I = 600, N = 399, 5(0.8 • 10" 3 ) « 1, see 11] for units. 



(fl"3*|l . not presented in Fig. [21 fits nicely to exact calcu- 
lation as long as t < t s . Then, the mean field approach 
fails to reproduce exact calculations at the front edge, 
and non-physical atom oscillations appear. We found 
that the QP term, playing a role when a spatial scale of 
density variations (s) becomes comparable to the healing 
length £, is responsible for their appearance. As t tends 
to t s we find s — > £ (Fig. [21 • Therefore, we find that the 
dynamical mean-field approach Q works nicely as long 
as s > ( and fails when s ~ £. In this way we bring 
a quantitative criteria into discussion 0, 0] concerning 
applicability of the Kolomeisky and Straley equation (0J 
to Tonks (Fermi) systems. 

Although, the derivative of density does not exactly 
go to infinity, the density profile becomes very steep at 
the front edge - Fig. [21 To look more quantitatively at 
shocks' dynamics we define the symmetry coefficient A as 

f Xm dx[p N (x) - Pn{1)]/ Jo™ dx[p N (x) - Pn(1)], i-e. as a 
ratio between the number of atoms in the front and rear 
impulse's parts with respect to the impulse maximum 
being at x m > 0. A is close to zero for a well devel- 
oped shock wave profile, and equals to one for symmetric 
perturbations. 

In the first stage of evolution, the profile undergoes 
sclf-steepening dynamics so that A decreases reaching the 
lowest value approximately at t = t s = 81.6 - Fig. [21[2t>- 
Then A(i) decreases a little and becomes flat - impulse 
propagates roughly without change of shape until t w 
150 (Fig. [21 Eb)- In the course of further evolution a 
shock wave front gradually disappears and the density 
profile becomes more and more symmetric (A — > 1) - Fig. 
[3i. Therefore, we conclude that there are three stages of 
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FIG. 3: Plot (a): further time evolution of the density profile 
from Fig. |2] Subsequent profiles correspond to t = 164, 262.4, 
344.4. Plot (b): symmetry coefficient A(t). Dotted line shows 
expected moment of the shock formation. Both plots present 
exact iV-body results 0- Parameters: as in Fig. H see H3 
for units. 
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FIG. 4: The exact many-body time evolution of the density 
profile of N = 399 fermions in a harmonic trap [§]. Sub- 
sequent profiles correspond to t = 0, 47r/120, 97t/120 (27r 
is a trap period) |T^|. Notice that for dark perturbations 
back edge self-steepens instead of the front one. Parameters: 
a = 0.75, u = 300. 



It is worthwhile to point out that the described split- 
ting of the initial perturbation has been already experi- 
mentally observed in a Bose-Einstein condensate (BEC) 
[To). An extension of our calculations to Bose gases pro- 
vides a theoretical description of a shock formation in a 
BEC. These results are presented in 16]. 

To summarize, we have studied in detail dynamics of 
density wave packets in the Fermi (Tonks) cloud. We 
have derived a simple analytical expression for the speed 
of propagation of Gaussian-like density wave packets hav- 
ing arbitrary amplitude. This result can be crucial for 
determination of the sound velocity from experimental 
data. We have also described analytically the process of 
a shock wave formation, and proposed a simple experi- 
mental setup for its observation. This prediction can be 
directly verified in a single run of experiment. Finally, we 
have found that the "quantum" shock wave front blows 
up at some point instead of becoming more and more 
steep in the course of free time evolution. The presence 
of shock waves in three-dimensional configurations can 
be an interesting subject for future studies. Indeed, as 
pointed out in jl3| , the propagation of density perturba- 
tions in a three-dimensional Fermi gas can differ signif- 
icantly from a quasi-one-dimensional case. The hydro- 
dynamical approach proposed in supported by exact 
many-body calculations can be used for such investiga- 
tions. 

I would like to acknowledge discussions with M. 
Brewczyk, Z. P. Karkuszewski, K. Sacha, and J. Za- 
krzewski. Work supported by KBN project 2 P03B 124 
22. 



evolution: the shock wave formation, propagation of a 
shock-like impulse, and the impulse explosion leading to 
broadening of the density profile. 

Finally, let us comment on the dynamics of density per- 
turbations in a harmonically trapped case: V(xi) = xf/2 
(Q. The division of an initial density perturbation into 
two similar pieces happens independently of the location 
of the laser beam focus. Just after the laser turn off 
the two perturbations travel in opposite directions to the 
edge of the cloud and come back to join after half-trap- 
period. At this moment the density profile is a mirror 
image of the initial shape. Then another splitting takes 
place and the dynamics repeats itself. As in the homo- 
geneous case, perturbations undergo self-steepening dy- 
namics leading to formation of a shock- like front, which 
then blows up - Fig. 0] 

Quasi-one-dimensional atomic traps form very promis- 
ing systems for observation of shock waves. Indeed, their 
reduced dimensionality avoids dispersion of the impulse 
amplitude due to multidimensional propagation. The pe- 
riodicity of shocks formation in a harmonic trap should 
allow for their multiple in situ observations [l(i| in a sin- 
gle run of an experiment. 
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